\chapter{Code Development}

As discussed previously, we intend to use the previous code from Martijn and work under the same framework. In order to advance one time step, the current process is shown in the following flowchart. Changing from Explicit to Implicit solver will require changes in the flowchart and subsequently solving different set of equations. We first start with the simple Test Case discussed in the previous Chapter, for which we have the benchmark solution available from MATLAB. Chorin's Projection scheme will be implemented in C++ for the Two-Dimensional Case. As we use the same framework, we should be able to plugin Martijn's RRB solver into our code without any modification. The combined solver in the step 1 of the Chorin's projection scheme will require a different strategy for solving, along with a discussion on different Preconditioners. This will be discussed in detail in following sections. 


\section{Combined Solver - GMRES/ BiCGStab}


For the first step of the Chorin's Projection Scheme, the system of equation to be solved is given by :

 \begin{equation}
  \dot{\textbf{q}} + L\textbf{q} =0,
 \end{equation}

with $\bq= \begin{bmatrix} \vec{\zeta} \\ \vec{\varphi} \end{bmatrix}$ and $\dot{\bq}$ its time derivative. 
The matrix $L = \begin{bmatrix} S_{\zeta \zeta} & S_{\zeta \varphi} \\ S_{\varphi \zeta} & S_{\varphi \varphi}\end{bmatrix}$ is the spatial discretization matrix as discussed in Chapter 2.

Implicit Trapezoidal method is used to integrate the system of ODE's in time given above. 

\begin{equation}
 (\dfrac{I}{\dt} + \beta L) q^{n+1} =  (\dfrac{I}{\dt} - (1-\beta) L) q^{n} 
\end{equation}


\subsection{Preconditioners}

Each GMRES or BiCGStab iteration requires the preconditioning step. The Combined preconditioner $M$ is given by:

$M = \begin{bmatrix} P_{\zeta \zeta} & Q_{\zeta \varphi} \\ Q_{\varphi \zeta} & P_{\varphi \varphi}\end{bmatrix}$ 

where $P$ represents the RRB preconditioner constructed over the block matrix $S$ , and $Q$ represents the off-diagonal blocks in the preconditioner.

Various preconditioners can be studies based on how $Q$ is represented.


$M_1 = \begin{bmatrix} P_{\zeta \zeta} & 0 \\ 0 & P_{\varphi \varphi}\end{bmatrix}$ is the simplest of the preconditioner. Solving the equation $Mz=r$ can be split into

\begin{equation}
 \begin{bmatrix} P_{\zeta \zeta} & 0\\ 0 & P_{\varphi \varphi}\end{bmatrix} \begin{bmatrix} z_{\zeta} \\ z_{\varphi} \end{bmatrix} = \begin{bmatrix} r_{\zeta} \\ r_{\varphi} \end{bmatrix} 
\end{equation}

which can be split into solving $P_{\zeta \zeta} z_{\zeta} = r_{\zeta}$ and $P_{\varphi \varphi} z_{\varphi} = r_{\varphi}$ :

\begin{equation}
 \begin{bmatrix} P_{\zeta \zeta} & 0\\ 0 & P_{\varphi \varphi}\end{bmatrix} \begin{bmatrix} z_{\zeta} \\ z_{\varphi} \end{bmatrix} = \begin{bmatrix} r_{\zeta} \\ r_{\varphi} \end{bmatrix} 
\end{equation}

$M_2 = \begin{bmatrix} P_{\zeta \zeta} & 0 \\ S_{\varphi \zeta} & P_{\varphi \varphi}\end{bmatrix}$ 

$M_3 = \begin{bmatrix} P_{\zeta \zeta} & S_{\zeta \varphi} \\ 0 & P_{\varphi \varphi}\end{bmatrix}$ 



\subsection{Preconditioner Setup}

\subsection{Preconditioning Solve}

\subsection{Impact of different Preconditioners on Spectrum}



